* PRoject: Ghana school feeding and child learning
* EA< 25 Sept 20202
* this dofile presents descriptive stats and other checks that were performed but only reported as text

******************** DESCRIPTIVE STATS CODES ******************	
	 
	 ***** TABLE 1: DESCRIPTIVE STATISTICS AND BALANCE COVARIATES AT BASELINE: TABLE PRESENTED RELATES TO THE FULL BASELINE SAMPLE
	use "$output/combined_data.dta", clear
	
	drop if child_panel ==0
	
	*Table balance of covariates
	
	reg agemo arm2 if wave==0 , vce(cluster locid)
	est sto h1
	
	reg male arm2 if wave==0 , vce(cluster locid)
	est sto h2
	
	reg enrol arm2 if wave==0 , vce(cluster locid)
	est sto h3
	
	reg grade_att arm2 if wave==0 , vce(cluster locid)
	est sto h4
	
	reg rep_grade arm2 if wave==0, vce(cluster locid)
	est sto h5
	
	reg absent arm2 if wave==0 , vce(cluster locid)
	est sto h6
	
	reg private arm2 if wave==0 , vce(cluster locid)
	est sto h7 
	
	reg haz arm2 if wave==0 , vce(cluster locid)
	est sto h8
	
	reg n_target arm2 if wave==0 , vce(cluster locid)
	est sto h9
	
	reg n_under5 arm2 if wave==0 , vce(cluster locid)
	est sto h10
	
	reg hhsize arm2 if wave==0 , vce(cluster locid)
	est sto h11
	
	reg headmale arm2 if wave==0 , vce(cluster locid)
	est sto h12
	
	reg mumage arm2 if wave==0 , vce(cluster locid)
	est sto h13
	
	reg mumedu arm2 if wave==0 , vce(cluster locid)
	est sto h14
	
	reg wi arm2 if wave==0 , vce(cluster locid)
	est sto h15
	
	reg sold arm2 if wave==0 , vce(cluster locid)
	est sto h16
	
	reg pcexp arm2 if wave==0 , vce(cluster locid)
	est sto h17
	
	reg livestock arm2 if wave==0 , vce(cluster locid)
	est sto h18
	
	reg urban arm2 if wave==0 , vce(cluster locid)
	est sto h19
	
	reg north arm2 if wave==0 , vce(cluster locid)
	est sto h20
	
	reg sick arm2 if wave==0 , vce(cluster locid)
	est sto h21
	
	reg headage arm2 if wave==0 , vce(cluster locid)
	est sto h21
	
	outreg2 [h*] using "$output/in_single3.xls", replace dec(3) stats(coef se) label 
	restore
	
	*** TABLE 2. REGRESSION OF CORRELATES OF REMAINING IN THE PANEL 
	
	est sto clear
	use "$output/combined_data.dta", clear
	
	duplicates drop qid wave, force
	***  Start by looking at the household attrition
	gen hh1= 1 if wave==0
	gen hh2=1 if wave ==1
	bys qid: egen hh_r1 =min(hh1)
	bys qid: egen hh_r2=min(hh2)
	bys qid: gen hh_panel = 0 if hh_r1==1 & hh_r2==. 
	bys qid: replace hh_panel =1 if hh_r1 ==1 & hh_r2==1
	
	***** attrition at the household level
	keep if wave==0
	collapse hh_panel arm2 locid n_target n_under5 hhsize headmale headage wi sold pcexp  mumage mumedu land_use livestock  north , by(qid)
	
	xi:logit hh_panel i.arm2 , vce(cluster locid) or
	est sto h1
	
	xi:logit hh_panel i.arm2 n_target n_under5 hhsize headmale headage wi sold pcexp livestock land_use north, vce(cluster locid) or
	est sto h2
	
	outreg2 [h*] using "$output/in_single3.xls", replace dec(3) stats(coef se) label eform
	
	***** attrition at the household level (check if differs by arm)
	keep if wave==0
	collapse hh_panel arm locid n_target n_under5 hhsize headmale headage wi sold pcexp  mumage mumedu land_use livestock  north , by(qid)
	
	xi:logit hh_panel i.arm , vce(cluster locid) or
	est sto h1
	
	xi:logit hh_panel i.arm n_target n_under5 hhsize headmale headage wi sold pcexp livestock land_use north, vce(cluster locid) or
	est sto h2
	
	outreg2 [h*] using "$output/in_single3.xls", replace dec(3) stats(coef se) label eform
	

	****** now look at child level attrition (proper table 2)
	** keep only panel children
	use "$output/combined_data.dta", clear
	***** new way to do attrition (as in Banerjee et al 2015, Science) ****** 7 MAY 2018
	
	est sto clear
	reg child_panel arm2 if wave==0, vce(cluster locid) 
	est sto h1
	
	reg child_panel arm2 zmaths zlit if wave==0, vce(cluster locid) 
	est sto h2
	
	reg child_panel arm2 male##arm2 north##arm2 poor_r1##arm2 if wave==0 , vce(cluster locid) 
	est sto h3
	
	outreg2 [h*] using "$output/in_single3.xls", replace dec(3) stats(coef se) label 
	
	
	** test if this is the same by three arms (noted in a footnote)
	
	use "$output/combined_data.dta", clear
	est sto clear
	xi:reg child_panel i.arm if wave==0, vce(cluster locid) 
	est sto h1
		
	xi:reg child_panel i.arm agemo male sick haz enrol n_target n_under5 hhsize headmale headage wi sold pcexp livestock land_use north if wave==0, vce(cluster locid) 
	est sto h2
	
	outreg2 [h*] using "$output/in_single3.xls", replace dec(3) stats(coef se) label 
	
	
	*** there is some non-differential attrition by some region and poor
	* investigate more about probability attrition poor in treatment
	
	ttest child_panel if wave==0 & poor_r1==1, by(arm2)
	
		*** TABLE 3. Descriptive statistics of outcomes at baseline and balance of outcomes at baseline tests +
	* descriptive statistics at endline (note: this is on longitudinal sample )
	use "$output/combined_data.dta", clear

	
	drop if child_panel ==0
	/*
	*********** SAMPLE AFTER ATTRITION *************
	tab arm2 wave

               |         wave
          arm2 |         0          1 |     Total
---------------+----------------------+----------
       Control |     1,502      1,502 |     3,004 
School feeding |     1,668      1,668 |     3,336 
---------------+----------------------+----------
         Total |     3,170      3,170 |     6,340 
*/
	
	***** NOTE: full longitudinal sample is 3,170, but for 187 children we dont have collected the values of cognitive scores
	estpost tabstat maths lit if wave==0 & arm2==0 , statistics(mean sd) columns(statistics) 
	
	esttab using "$output/IN_estimates3.csv", replace cells(mean(fmt(2)) sd(par)) nostar unstack obslast ///
	nonote nomtitle nonumber label title(YC - Round 2) 
	
	estpost tabstat maths lit if wave==0 & arm2==1 , statistics(mean sd) columns(statistics) 
	
	esttab using "$output/IN_estimates3.csv", replace cells(mean(fmt(2)) sd(par)) nostar unstack obslast ///
	nonote nomtitle nonumber label title(YC - Round 2) 
	
	reg maths arm2 agemo if wave==0 , vce(cluster locid)
	est sto h1
	
	reg lit arm2 agemo if wave==0 , vce(cluster locid)
	est sto h2
	
	outreg2 [h*] using "$output/in_single3.xls", replace dec(3) stats(coef se) label 
	
	* to see whether they become sign
	
	reg zmaths arm2 agemo if wave==1 , vce(cluster locid)
	est sto h1
	
	reg zlit arm2 agemo if wave==1 , vce(cluster locid)
	est sto h2
	
	outreg2 [h*] using "$output/in_single3.xls", replace dec(3) stats(coef se) label 
	
	*** Figure 3: raw scores by gender, household poverty and northern (this is for full long sample)
	use "$output/combined_data.dta", clear
	drop if child_panel ==0
	*** full sample maths and lit
	
	kdensity math if wave == 0 & arm2==0 , lwidth(medthick) lcolor(black) addplot(kdensity math if wave == 0 & arm2==1, lwidth(medthick) lpattern(dash)|| ///
		kdensity math if wave == 1 & arm2==0 , lwidth(thick) lcolor(black) lpattern(dot) || kdensity math if wave == 1 & arm2==1 , lwidth(medthick) lpattern(longdash_dot)) ///
		legend(ring(0) pos(2) label(1 "Control baseline") label(2 "School feeding baseline")  ///
				label(3 "Control endline") label(4 "School feeding endline") ) graphregion(color(white)) title("") scheme(s2mono) ///
				xtitle("Math raw score")

	 graph save Graph "$output/graphsoct2019/maths_all_sample.gph", replace

	kdensity lit if wave == 0 & arm2==0 , lwidth(medthick) lcolor(black) addplot(kdensity lit if wave == 0 & arm2==1, lwidth(medthick) lpattern(dash) || ///
		kdensity lit if wave == 1 & arm2==0 , lwidth(thick) lcolor(black) lpattern(dot) || kdensity lit if wave == 1 & arm2==1 , lwidth(medthick) lpattern(longdash_dot)) ///
		legend(ring(0) pos(2) label(1 "Control baseline") label(2 "School feeding baseline")  ///
				label(3 "Control endline") label(4 "School feeding endline") ) graphregion(color(white)) title("") scheme(s2mono)

	graph save Graph "$output/graphsoct2019/lit_all_sample.gph", replace
		
	grc1leg "$output/graphsoct2019/maths_all_sample.gph"  "$output/graphsoct2019/lit_all_sample.gph" , legendfrom("$output/graphsoct2019/lit_all_sample.gph") graphregion(color(white)) 
	graph save Graph "$output/graphsoct2019/mathlit_all_sample.gph", replace

	**** with z scores
	
	kdensity zmaths if wave == 0 & arm2==0 , lwidth(medthick) lcolor(black) addplot(kdensity zmaths if wave == 0 & arm2==1, lwidth(medthick) lpattern(dash) || ///
		kdensity zmaths if wave == 1 & arm2==0 , lwidth(thick) lcolor(black) lpattern(dot) || kdensity zmaths if wave == 1 & arm2==1 , lwidth(medthick) lpattern(longdash_dot)) ///
		legend(ring(0) pos(2) label(1 "Control baseline") label(2 "School feeding baseline")  ///
				label(3 "Control endline") label(4 "School feeding endline") ) graphregion(color(white)) title("") scheme(s2mono) ///
				xtitle("Age-standardized math score")

	 graph save Graph "$output/graphsoct2019/zmaths_all_sample.gph", replace

	kdensity zlit if wave == 0 & arm2==0 , lwidth(medthick) lcolor(black) addplot(kdensity zlit if wave == 0 & arm2==1, lwidth(medthick)  lpattern(dash) || ///
		kdensity zlit if wave == 1 & arm2==0 , lwidth(thick) lcolor(black) lpattern(dot) || kdensity zlit if wave == 1 & arm2==1 , lwidth(medthick) lpattern(longdash_dot)) ///
		legend(ring(0) pos(2) label(1 "Control baseline") label(2 "School feeding baseline")  ///
				label(3 "Control endline") label(4 "School feeding endline") ) graphregion(color(white)) title("") scheme(s2mono) ///
				xtitle("Age-standardized literacy score")

	graph save Graph "$output/graphsoct2019/zlit_all_sample.gph", replace
		
	grc1leg "$output/graphsoct2019/zmaths_all_sample.gph"  "$output/graphsoct2019/zlit_all_sample.gph" , legendfrom("$output/graphsoct2019/zlit_all_sample.gph") graphregion(color(white)) 
	graph save Graph "$output/graphsoct2019/zmathlit_all_sample.gph", replace

	******************************* FIGURE 4 **********************************
	*************+ BY GENDER		
	
	**** Maths and lit by gender 
	kdensity zmaths if wave == 1 & arm2==0 & male==1 , lwidth(medthick) lcolor(black) addplot(kdensity zmaths if wave == 1 & arm2==0 & male==0, lwidth(medthick) lpattern(dash) || ///
		kdensity zmaths if wave == 1 & arm2==1 & male==1, lwidth(thick) lcolor(black) lpattern(dot) || kdensity zmaths if wave == 1 & arm2==1 & male==0 , lwidth(medthick) lpattern(longdash_dot)) ///
		legend(ring(0) pos(2) label(1 "Control boys") label(2 "Control girls")  ///
				label(3 "School feeding boys") label(4 "School feeding girls") ) graphregion(color(white)) title("") ///
				xtitle("Age-standardized math score") scheme(s2mono)
	
	 graph save Graph "$output/graphsoct2019/zmaths_endline_gender_all.gph", replace

	kdensity zlit if wave == 1 & arm2==0 & male==1 , lwidth(medthick) lcolor(black) addplot(kdensity zlit if wave == 1 & arm2==0 & male==0, lwidth(medthick) lpattern(dash) || ///
		kdensity zlit if wave == 1 & arm2==1 & male==1, lwidth(thick) lcolor(black) lpattern(dot) || kdensity zlit if wave == 1 & arm2==1 & male==0 , lwidth(medthick) lpattern(longdash_dot)) ///
		legend(ring(0) pos(2) label(1 "Control boys") label(2 "Control girls")  ///
				label(3 "School feeding boys") label(4 "School feeding girls") ) graphregion(color(white)) title("") ///
				xtitle("Age-standardized literacy score") scheme(s2mono)
	
	 graph save Graph "$output/graphsoct2019/zlit_endline_gender_all.gph", replace
	 
	 grc1leg "$output/graphsoct2019/zmaths_endline_gender_all.gph" "$output/graphsoct2019/zlit_endline_gender_all.gph" , legendfrom("$output/graphsoct2019/zlit_endline_gender_all.gph") graphregion(color(white)) 
	 graph save Graph "$output/graphsoct2019/zmaths_lit_gender.gph", replace
	 
	 ***** BY POVERTY STATUS ***********
	
	kdensity zmaths if wave == 1 & arm2==0 & poor_r1==0 , lwidth(medthick) lcolor(black) addplot(kdensity zmaths if wave == 1 & arm2==0 & poor_r1==1, lwidth(medthick) lpattern(dash) || ///
		kdensity zmaths if wave == 1 & arm2==1 & poor_r1==0, lwidth(thick) lcolor(black) lpattern(dot) || kdensity zmaths if wave == 1 & arm2==1 & poor_r1==1 , lwidth(medthick) lpattern(longdash_dot)) ///
		legend(ring(0) pos(2) label(1 "Control - Above poverty line") label(2 "Control - below poverty line")  ///
				label(3 "School feeding - above poverty line") label(4 "School feeding - below poverty line") ) graphregion(color(white)) title("") ///
				xtitle("Age-standardized math score") scheme(s2mono)
	
	 graph save Graph "$output/graphsoct2019/zmaths_endline_poverty_all.gph", replace
	
	kdensity zlit if wave == 1 & arm2==0 & poor_r1==0 , lwidth(medthick) lcolor(black) addplot(kdensity zlit if wave == 1 & arm2==0 & poor_r1==1, lwidth(medthick) lpattern(dash) || ///
		kdensity zlit if wave == 1 & arm2==1 & poor_r1==0, lwidth(thick) lcolor(black) lpattern(dot) || kdensity zlit if wave == 1 & arm2==1 & poor_r1==1 , lwidth(medthick) lpattern(longdash_dot)) ///
		legend(ring(0) pos(2) label(1 "Control - Above poverty line") label(2 "Control - below poverty line")  ///
				label(3 "School feeding - above poverty line") label(4 "School feeding - below poverty line") ) graphregion(color(white)) title("") ///
				xtitle("Age-standardized math score") scheme(s2mono)
				
	
	 graph save Graph "$output/graphsoct2019/zlit_endline_poverty_all.gph", replace
	
	grc1leg "$output/graphsoct2019/zmaths_endline_poverty_all.gph" "$output/graphsoct2019/zlit_endline_poverty_all.gph", legendfrom("$output/graphsoct2019/zlit_endline_poverty_all.gph") graphregion(color(white)) 
	 graph save Graph "$output/graphsoct2019/zmathlit_poverty.gph", replace
	 
	**************** BY REGION ***************
	kdensity zmaths if wave == 1 & arm2==0 & north==0 , lwidth(medthick) lcolor(black) addplot(kdensity zmaths if wave == 1 & arm2==0 & north==1, lwidth(medthick) lpattern(dash) || ///
		kdensity zmaths if wave == 1 & arm2==1 & north==0, lwidth(thick) lcolor(black) lpattern(dot) || kdensity zmaths if wave == 1 & arm2==1 & north==1 , lwidth(medthick) lpattern(longdash_dot)) ///
		legend(ring(0) pos(2) label(1 "Control - Southern regions") label(2 "Control - Northern regions")  ///
				label(3 "School feeding - Southern regions") label(4 "School feeding - Northern regions") ) graphregion(color(white)) title("") ///
			xtitle("Age-standardized math score") scheme(s2mono)	
	
	 graph save Graph "$output/graphsoct2019/maths_endline_north_all.gph", replace
	
	kdensity zlit if wave == 1 & arm2==0 & north==0 , lwidth(medthick) lcolor(black) addplot(kdensity zlit if wave == 1 & arm2==0 & north==1, lwidth(medthick) lpattern(dash)  || ///
		kdensity zlit if wave == 1 & arm2==1 & north==0, lwidth(thick) lcolor(black) lpattern(dot) || kdensity zlit if wave == 1 & arm2==1 & north==1 , lwidth(medthick) lpattern(longdash_dot)) ///
		legend(ring(0) pos(2) label(1 "Control - Southern regions") label(2 "Control - Northern regions")  ///
				label(3 "School feeding - Southern regions") label(4 "School feeding - Northern regions") ) graphregion(color(white)) title("") ///
				xtitle("Age-standardized math score") scheme(s2mono)
	
	graph save Graph "$output/graphsoct2019/lit_endline_north_all.gph", replace
	
	grc1leg "$output/graphsoct2019/maths_endline_north_all.gph" "$output/graphsoct2019/lit_endline_north_all.gph", legendfrom("$output/graphsoct2019/lit_endline_north_all.gph") graphregion(color(white)) 
	 graph save Graph "$output/graphsoct2019/mathlit_north.gph", replace
	 	
	* Correlation test scores in other datasets 
	* i use young lives PPVT
	
	use  "/Users/betta/Dropbox/Sackler_PRPR special issue/output/dataset_pnrp.dta", clear
	keep childid round ppvt_raw   
	reshape wide ppvt_raw, i(childid) j(round)
	pwcorr ppvt_raw2 ppvt_raw3 ppvt_raw4 ppvt_raw5, star(.001)
	
	 ***** appendix: correlates of compliance/uptake
	use "$output/combined_data.dta", clear
	drop if child_panel ==0
 
est sto clear
	xi:logit sf agemo male i.exp_c poor_r1 i.north  if wave ==1 & arm2==1, vce(cluster locid) or
	est sto h1
	
	xi:logit sf agemo male i.exp_c poor_r1 zmaths_r1 zlit_r1 zraven_r1 zdigit_r1 grade_r1  i.north  if wave ==1 & arm2==1, vce(cluster locid) or
	est sto h2
	
	xi:logit sf agemo male i.exp_c poor_r1 i.north zmaths_r1 zlit_r1 zraven_r1 zdigit_r1 grade_r1 i.arm  if wave ==1 & arm2==1, vce(cluster locid) or
	est sto h3
	
	xi:logit sf agemo male i.exp_c poor_r1 zmaths_r1 zlit_r1 zraven_r1 zdigit_r1 grade_r1 i.north primary n_child_sf private haz_r1 if wave ==1 & arm2==1, vce(cluster locid) or
	est sto h4
	
	outreg2 [h*] using "$output/in_single3.xls", replace dec(3) stats(coef se) label eform
	
	
	* check if proportion of children going to private changed between baseline and endline
	use "$output/combined_data.dta", clear
	
	
	drop if child_panel ==0
	
	keep u_id child_panel zmaths zlit zraven zdigit learn_ind cog_ind all_index arm arm2 wave region agemo locid male north poor_r1 ///
	enrol primary days_attend grade rep_grade private sch_time stu_time hrs_care hrs_work  hrs_lei exp_c arm sick haemo ///
	not_working day_hhwork day_other_biz day_fam_biz day_work_tot num_meal n_break nobrk dd9 
	
	reshape wide child_panel agemo zmaths zlit zraven zdigit learn_ind cog_ind all_index arm2 region locid male north poor_r1 ///
	enrol primary days_attend grade rep_grade private sch_time stu_time hrs_care hrs_work  hrs_lei exp_c arm sick haemo ///
	not_working day_hhwork day_other_biz day_fam_biz day_work_tot num_meal n_break nobrk dd9 , i(u_id) j(wave)
	ssss
	
	
	reg agemo arm2 if wave==0 & poor_r1==1 , vce(cluster locid)
	
	* check if proportion of private changed between baseline and endline (e.g. switching from private to public)

	xi: reg private1 private0##arm20 i.region0 , cluster(locid0) 
	
	*** data on uptake reported in the paper
	
	su sf if wave==1 & arm2==1 
	*61%
	su sf if wave ==1 & arm2 ==1 & private ==0 & primary==1
	* 82.2
	
	* gen mean uptake by community
	bys locid: egen mean_uptake = mean(sf) if wave ==1 & arm2 ==1 & private ==0 & primary==1
	su mean_uptake
	histogram mean_uptake
	
	. tab mean_uptake

	/*
mean_uptake	Freq.	Percent	Cum.
			
0	25	2.63	2.63
.0909091	11	1.16	3.79
.1	10	1.05	4.84
.2666667	15	1.58	6.41
.3333333	18	1.89	8.31
.4444444	9	0.95	9.25
.4615385	13	1.37	10.62
.4827586	29	3.05	13.67
.5	10	1.05	14.72
.5263158	19	2.00	16.72
.5714286	14	1.47	18.19
.7307692	26	2.73	20.93
.75	12	1.26	22.19
.7777778	9	0.95	23.13
.7843137	51	5.36	28.50
.7857143	28	2.94	31.44
.7894737	19	2.00	33.44
.8461538	26	2.73	36.17
.85	20	2.10	38.28
.8888889	18	1.89	40.17
.9047619	21	2.21	42.38
.90625	32	3.36	45.74
.9166667	24	2.52	48.26
.9333333	15	1.58	49.84
.9375	16	1.68	51.52
.9411765	17	1.79	53.31
.9615385	26	2.73	56.05
.9642857	28	2.94	58.99
.96875	32	3.36	62.36
1	358	37.64	100.00
			
Total	951	100.00
*/

	************ ADDITIONAL INFO REPORTED IN THE PAPER **************
	
	use "$output/combined_data.dta", clear
	
	drop if child_panel ==0
	
	keep u_id child_panel zmaths zlit zraven zdigit learn_ind cog_ind all_index arm arm2 wave region agemo locid male north poor_r1 ///
	enrol primary days_attend grade rep_grade private sch_time stu_time hrs_care hrs_work  hrs_lei exp_c arm sick haemo ///
	not_working day_hhwork day_other_biz day_fam_biz day_work_tot num_meal n_break nobrk dd9 
	
	reshape wide child_panel agemo zmaths zlit zraven zdigit learn_ind cog_ind all_index arm2 region locid male north poor_r1 ///
	enrol primary days_attend grade rep_grade private sch_time stu_time hrs_care hrs_work  hrs_lei exp_c arm sick haemo ///
	not_working day_hhwork day_other_biz day_fam_biz day_work_tot num_meal n_break nobrk dd9 , i(u_id) j(wave)
	ssss
	
	****  check correlations of variables (reported in the paper)

	pwcorr zmaths1 zmaths0, star(0.01)
	pwcorr zlit1 zlit0, star(0.01)
	pwcorr zdigit1 zdigit0, star(0.01)
	pwcorr zraven1 zraven0, star(0.01)
	pwcorr learn_ind1 learn_ind0, star(0.01)
	pwcorr cog_ind1 cog_ind0, star(0.01)
	pwcorr all_index1 all_index0 , star(0.01)
	
	pwcorr enrol1 enrol0 , star(0.01)
	pwcorr primary1 primary0, star(0.01)
	pwcorr days_attend1 days_attend0, star(0.01)
	pwcorr grade1 grade0 , star(0.01)
	pwcorr rep_grade1 rep_grade0, star(0.01)

	* check if missing depends on treatment assignment *

	* gen variable no missing 

	gen no_missing1 =1 if zmaths1!=. & zlit1!=. & zdigit1!=. & zraven1!=. 
	gen no_milear=1 if zmaths1!=. & zlit1!=.
	gen no_micog=1 if zdigit1!=. & zraven1!=. 
	
	gen missing_math =1 if zmaths1==.
		replace missing_math=0 if zmaths1!=.
		
	foreach v of varlist zlit1 zdigit1 zraven1 {
		gen missing_`v'= 1 if `v'==.
			replace missing_`v' =0 if `v'!=.
			}
	
	est sto clear
	reg missing_math arm20 , cluster(locid0)
	est sto h1
	
	reg missing_zlit1 arm20 , cluster(locid0)
	est sto h2
	
	reg missing_zdigit1 arm20 , cluster(locid0)
	est sto h3
	
	reg missing_zraven1 arm20 , cluster(locid0)
	est sto h4
	
	outreg2 [h*] using "$output/in_single3.xls", replace dec(3) stats(coef se) label sortvar(arm2) 
	
	

	
	
	
